#!/usr/bin/R
	
#set model parameters
	t.horizon=100 #simulation length (hrs)
	dt=.1 #timestep length (hrs)
	ET=0 #evapotranspiration rate (m/hr) ###consider functional relationship to temperature???###	
	PP=0 #precipitation rate (m/hr) ###forecasted???###
	grav=32.2 #ft/s^2

#intialize variables
	source('pond.params.4.R') #must source this file BEFORE variable initialization
	n.ponds=nrow(pond)
	n.boxs=nrow(box)
	n.times=as.integer(t.horizon/dt) # number of timesteps
	h=array(0,dim=c(n.ponds,n.times)) #water surface elev above MSL
	plant.inflow=array(0,dim=n.times) #inflow to first pond (cms)
	pond.inflow=array(0,dim=c(n.ponds,n.times)) #net inflow to each pond (cms)
	pond.outflow=array(0,dim=c(n.ponds,n.times)) #net outflow from each pond (cms)
	#q=array(0,dim=c(n.boxs,n.times)) #flow past each connection (cms)
	#q.temp=array(0,dim=n.boxs) #dummy inflow for RK routine
	K=array(0,dim=c(3,n.ponds)) #dummy variables for RK routine
	h[,1]=as.numeric(pond[,p.surf]) #set initial conditions
	area=as.numeric(pond[,p.area]) #set pond areas
	#colors##########################################
	library(fields,verbose=FALSE)
        my.colors=array(tim.colors(),dim=(n.times+1))
	#colors##########################################
	source('pondSim.module.3.R') #must source this file AFTER variable initialization	
		
#perform simulation
	#png(file='images/1.png',width=500,height=400,units='px',bg='transparent')
	barplot(h[,1],col=my.colors[1],ylim=c(0,15)) #plot initial conditions
	#dev.off()
	if (n.times==1) stop()
	for (t in 1:(n.times-1)) {
		#RK STEP 1
		foo=find.flow(h[,t])
			pond.inflow=foo[[1]]+(PP*area)
			pond.outflow=foo[[2]]+(ET*area)
			pond.inflow[1]=pond.inflow[1]+plant.inflow[t] #add plant inflow to first pond
			K[1,]=dt/area*(pond.inflow-pond.outflow)
			h.temp=h[,t]+(1/3)*K[1,]
			h.temp[h.temp<0]=0
		#RK STEP 2
		foo=find.flow(h.temp)
                        pond.inflow=foo[[1]]+(PP*area)
                        pond.outflow=foo[[2]]+(ET*area)
                        pond.inflow[1]=pond.inflow[1]+plant.inflow[t] #add plant inflow to first pond
                        K[2,]=dt/area*(pond.inflow-pond.outflow)
                        h.temp=h[,t]+(2/3)*K[2,]
			h.temp[h.temp<0]=0
		#RK STEP 3
		foo=find.flow(h.temp)
                        pond.inflow=foo[[1]]+(PP*area)
                        pond.outflow=foo[[2]]+(ET*area)
                        pond.inflow[1]=pond.inflow[1]+plant.inflow[t] #add plant inflow to first pond
                        K[3,]=dt/area*(pond.inflow-pond.outflow)
                        h.temp=h[,t]+(1/4)*K[1,]+(3/4)*K[3,]
			h.temp[h.temp<0]=0
		#RK STEP 4
		h[,t+1]=h.temp
		#png(file=paste('images/',t+1,'.png',setp=''),width=500,height=400,units='px',bg='transparent')
		barplot(h[,t+1],col=my.colors[t+1],ylim=c(0,15))
		#dev.off()
	}
